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Abstract 

n 

We consider a chain of Josepshon-junction rhombi (proposed originally in in quantum regime, 
and in the realistic case when charging effects are determined by junction capacitances. In the 
maximally frustrated case when magnetic flux through each rhombi <3? r is equal to one half of su- 
perconductive flux quantum <&q, Josepshon current is due to correlated transport of pairs of Cooper 
pairs, i.e. charge is quantized in units of 4e. Sufficiently strong deviation 5<& = |<3? r — <I>o / 2 1 > S& c 
from the maximally frustrated point brings the system back to usual 2e-quantized supercurrent. 
We present detailed analysis of Josepshon current in the fluctuation-dominated regime (sufficiently 
long chains) as function of the chain length, Ej/Eq ratio and flux deviation 5$. We provide 
estimates for the set of parameters optimized for the observation of 4e-super cur rent. 



I. INTRODUCTION 



Pairing of Cooperpairs in frustrated Josephson junction arrays was theoretically pro- 
posed recently [110,0] in the search of topologically protected nontrivial quantum liquid 
states. The simplest system where such a phenomenon could be observed was proposed by 
Doucot and Vidal in It consists of a chain of rhombi (each of them being small ring of 
4 superconductive islands connected by 4 Josephson junctions) placed into transverse mag- 
netic field, cf. Fig. 1. It was shown in jl| that in the fully frustrated case (i.e. magnetic 
flux through each rhombus $ = = t~ ) usual tunnelling of Cooper pairs along the chain 
is blocked due to destructive inteference of tunneling going through two paths within the 
same rhombus, while correlated 2-Cooper-pair transport survives. Evidently, experimental 
observation of such a phenomenon (detected as anomalous period |$o of the global super- 
current along the chain) would be very desirable. However, theoretical results of Ref. Q] refer 
to the situation when Coulomb energy is determined by self-capacitances Cq of individual 
superconductive islands, whereas in real submicron Josephson-j unction arrays capacitances 
of junctions C dominate, cf. e.g. j^J. In this paper we reconsider the model of Ref. 0] 
for the experimentally relevant situation C ^> Cq. This case is also simpler for theoretical 
treatment, since Lagragian of the system becomes a sum of terms, such that each of them 
belongs to individual rhombus only. The only source of coupling between different rhombi is 
the periodic boundary condition along the chain. The method to treat similar problem was 
developed recently by Matveev, Larkin and Glazman 0] (MLG). They considered simple 
chain of N Josephson junctions in the closed- ring geometry, and reduced calculation of su- 
percurrent in large- N limit to the solution of a Schroedinger equation for a particle moving 
in a periodic potential ~ cosx, with appropriate boundary condition. MLG assumed (we 
will do the same) that Josephson energy Ej of junctions is large compared to their charging 
energy Eq = e 2 /2C. We will generalize the MLG method in order to use it for the case 
of ring of rhombi. It will be shown that in our case fictitios particle of the MLG theory is 
still moving in the cos-like potential, but it acquires now large spin S = |iV, where N is 
the number of rhombi in the ring. In the maximally frustrated case |<3> r — $o/2| = 5$ = 
the ^-projection of the spin is an integral of motion, which should be chosen to minimize 
the total energy. As a result, S x = ±|iV and the whole problem reduces to the one studied 
by MLG up to trivial redefinition of parameters. In this situation ground-state energy and 
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supercurrent (which is proportional to derivative of the ground-state energy over total flux 
$ c ) are periodic function of $ c with period $o/2, i.e. 4e-transport takes place. Nonzero flux 
deviation produces longitudinal field h z coupled to the z-component of spin of fictitios 
particle, which acquires now nontrivial dynamics. We show that in the limit of sufficiently 
long rhombi chain the whole problem can be analyzed in terms of semiclassical dynamics of 
a particle with a large spin under spin-dependent potential barrier. In general, there are two 
tunnelling trajectories, one of them corresponds to usual 2e transport, whereas another to 
4e transport. Comparing actions of these trajectories for different we find critical flux 
deflection <5$ c as function of the ratio Ej/Ec 3> 1. 

The rest of the paper is organized as follows: in Sec. II we define our model and identify 
its classical states; in Sec. Ill we derive effective Hamiltonian which governs quantum phase 
slip processes and calculate supercurrent as function of the flux deflection <5$ c ; in Sec. IV we 
consider current-bias chain with I > I c and calculate voltage V(I) via the rate of incoherent 
quantum phase slips. Our conclusions and suggestions for the experiment are presented in 
Sec. V. Finally, in Appendix A somewhat tedious calculation of current-phase relation is 
presented. 

II. THE MODEL AND ITS CLASSICAL STATES 




N 1 2 

Figure 1: The chain of rhombi (IS i\ closed ring. 

We study a chain of N rhombi shown in Fig. ^ Each rhombi consists of four supercon- 
ductive islands connected by tunnel junctions with Josephson coupling energy Ej = hl®/2e; 
charging energy Eq is determined by capacitance C of junctions, Ec = e 2 /2C (we neglect 
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self-capacitances of islands which are assumed to be much smaller than C). Below we con- 
sider Josephson current along the chain of N ^> 1 rhombi and assume that the chain is of 
the ring shape, with total magnetic flux $ c inside the ring. We also denote by $ r the flux 
per elementary rhombus plaquette and define phases 7 and ip: 

where $ — h/2e is the superconducting flux quantum. We study the situation than <3> is 
close to $o/2, i.e. S — ip — n <C 1. 

Assuming that the charge transport through the system at 5 = is carried out by charges 
4e, we expect that in this case dependence of the current in the chain on the external flux 
$ c is periodic with period $o/2. Below we calculate the $o/2-periodic current at 5 = 0. 
We also show that at small 5 the current through the system has two components J 4e and 
I 2e with periods $o/2 and <E> respectively. The first component corresponds to the current 
carried by pairs of Cooper pairs and the second one corresponds to the single Cooper pair 
transport. At very small 5 the current dominates over J 2e . We will refer to this regime 
as to 4e-regime. At large enough S the opposite situation (2e-regime) is realized. We will 
find below the crossover point S c between these two regimes. 

The system is described by the following imaginary-time action 



1 (de^ : 
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Here the variable 6^ is the phase difference across the m-th junction in the n-th rhombus 
(see Fig. 1). Taking into account that each rhombus is pierced by flux $ r and the flux 
through the whole chain is $ c we derive the following additional conditions 

4 

J2e^ = V , n=l,2,..,iV, (3) 

m=l 

N 

EM 3) -a=7- (4) 



n=l 



In this paper we consider the case of strong coupling between grains Ej ^> Eq- This 
enables us to use semi-classical approximation for calculating the energy spectrum of the 
system. At E c = the phases 6*1"^ become classical variables and the energy states of the 
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chain can be found by minimizing the sum of Josephson energy terms in action (J2J). Let us 
introduce variables 8 n = —9n — @n \ where 9 n is the phase difference along the diagonal 
of the n-th rhombus. It is convenient to make minimization in two steps. First of all the 
Josephson energy of a single rhombus under the fixed flux through the rhombus and under 
the fixed phase difference 9 n is minimized. For the Josephson energy of the chain we then 
get for 5 <C 1: 



E = -2V2Ej + \ 6 <) cos (y - 



(5) 



sin/3 n = (l - |<) , oosp n = ±-L + 5 - a ^j . ( 6 ) 

Plus and minus signs in (jUJ) correspond to positive (resp. negative) values of cos|. Here 
we have introduced an important notation a z = sign sin# n . It can be easily shown that at 
5 = each individual rhombus has two classical ground states with equal energies. This 
states differ only in the sign of the superconducting current circulating around the plaquette 
which corresponds to the binary variable a z n . 

Now we have to minimize the energy (jSJ) with respect to phases On subject to the 
constrains (JIJ. Assuming 5 to be small and N to be large we get 

E m ,a « ^r(7 - kN/2 - tiS z - 2nm) 2 - V25S z Ej + Const. (7) 

Here m is an arbitrary integer (which has the same meaning as in the MLG paper [5]) and 

1 z oz * ~ N( P TIN N6 

L 7 = 7 + ^ = 7 + ^ + ^- (8) 
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In the above equation s z n can be considered as ^-projection of the "spin" ^ which describes 
binary degeneracy of states of the n-th rhombi. Then S z corresponds to the ^-projection of 
the total large spin S describing the whole rhombi chain. For clarity everywhere in this paper 
we will refer to the case of even number of rhombi. Then total spin S and eigenvalues of its 
projection S z are integer. Classical states of the chain are characterized by individual spin 
projections a z n for each rhombus, and by collective integer-valued variable m. We will denote 
these states by \m, {c^}) or \m,a >. Physically, classical state of the chain is characterized 
by the global current / along the chain, and by the signs of local currents flowing in each 
of N rhombi. Nonzero charging energy Eq provides quantum phase slips in each of 4iV 
Josephson junction; these processes mix different classical states leading to formation of the 
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Figure 2: The numerical factor A; as a function of Ej/E c . 

fully quantum ground state. Below we derive effective Hamiltonian acting on the space of 
classical states, and find ground-state energy E (j) and corresponding supercurrent. 

III. QUANTUM FLUCTUATIONS OF RHOMBI AND SUPERCURRENT 

We turn to analysis of quantum fluctuations of 6^ at finite E C - The most important 
type of these fluctuations involve an instanton (quantum phase slip, QPS), i.e. a trajectory 
that begin at one of the minima (|IJ) of the potential energy in action (J2J) at t = — oo and 
ends at another minima at t = +00. There are two kinds of trajectories: the first starts 
at \m, {&n}) and ends at |m, \p z n + 2<5 n fc}) for arbitrary 1 < k < N, o z k = —1, whereas the 
second start at \m, {cr^}) and end at \m + 1, \o z n — 2<5 nfc }) for arbitrary 1 < k < N, cr| = 1. 
Any trajectory of the first kind corresponds to QPS in 6% or whereas trajectory of 
the second kind corresponds to QPS 9^ or 6f\ Note that at 5 = and 7 = tv/2 all these 
trajectories starting at |m, {u^}) with 2m + S z = connect the minima with equal energies. 
Thus they are important for restoring symmetry of the system which is classically broken. 
Let us denote as v the amplitude of a QPS in one contact. At large N 3> 1 this amplitude 
does not differ from the "spin flip" amplitude for a single rhombus at $ r ~ $o/2. In this 
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approximation we can use result from Ref 
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k{E)E c ) l/ \^[-lM x j^-). (9) 



where k is a numerical factor of the order of unity. Comparison with direct numerical diag- 
onalization [6.] of of low-lying spectrum of a single frustrated rhombi shows that coefficient 
k grows from approximately 1.3 to 1.44 as the ratio Ej/Eq varies from 10 to infinity, cf. 
Fig. El The instantons account for the possibility of the system tunneling between different 
minima (JZJ) of the potential energy in action (J2J). The effect of the instantons on the ground 
state energy can be represented by a tight-binding Hamiltonian defined as 

#KK}) = £m,<xK{<}> + 

N 



fc=i 

N 

2v h+M^---^i>-4>---4+i><}) + 

k=l , o z k =\ 

N 

2v \ m ~ 1 '{ (T i>"-- cr fc-i'-^>--- or fc+i>^}> • ( 10 ) 

fc=l,<7*=-l 

To find the ground state energy £"(7) it is convenient to make Fourier transformation 
over variable m according to 

\x,a) = Y.^{ 2i ( 2m - 1 + S z + ^)x\\m,a), S 2 = ^<. (11) 

m ^ ^ '-^ n=l 

Note that not all vectors of our new basis (fTT|) are independent. It is easy to see from (fTT]) 
that \x + 7r/2, a) = e - t7 y+ mSZ + mN / 2 \ x , a). Considering any system's state of the form 
= Yl x a a ) \ x i a ) we should impose on the wavefunction ip(x, a) a twisted boundary 
condition 

Here we have introduced operator S z acting on the spin variables of the system according 
to standart rules. 

The resulting Schrodinger equation acquires the form 

d 2 i) 



dx 2 



+ (b - 2w cos 2x ■ S x + 2hS z )ip = 0, (13) 



where 

1QNE 6ANv , 8N5 , , 

b=—= , w = —= , h — — — . 14 

Note that symmetry group of the Hamiltonian corresponding to the equation (j!3|) includes 
transformations U n = e lwn ( s +N ^ 2 ^T nn /2- Here T a is operator of the translation over distance 
a along the x-axis. Equation (fE?j) shows that the parameter 7 specifies different irreducible 
representations of the symmetry group. 

The equations ()13ll2j) allow comprehensive analytical investigation for the case when 
the flux per single rhombus equals $o/2- For such a system h = and the Hamiltonian 
commutes with S x . However, variables x and S are still coupled due to boundary condition 
|T2|) . Therefore we look for wavefunction in the form 

il>(x, a) = \S X ) 4>{x) + e indz+i7lN/2 \S X ) <j>(x + tt/2) (15) 

Then for <fi{x) we get standard Mathieu equation: 

+ (b- 2gcos2x)0 = 0, (16) 



dx 2 

where q = wS x . Boundary condition (|T2*|) now reads as follows: 

(f)(x + tt) = e 2l ^<f){x) . (17) 

The ground state of the system defined by Eqs. IT?]) corresponds to maximal absolute 
value of S x , equal to N/2. In other words, in the ground state all rhombi in the chain 
are either in symmetric superpositions or in antisymmetric superpositions, of their double- 
degenerate classical states. Thus there are two degenerate eigenstates 

|0±) = Ux)\S x = N/2) and |0~) = <fr{x + n/2)\S x = -N/2) = Ux\0±) , (18) 

with the same lowest energy Eq. This degeneracy is a direct consequence of the fact that 
for h = (fully frustrated chain) the Hamiltonian has two symmetry operators S x and U\ 
which do not commute with each other (thus double-degeneracy refers to all states, not only 
to the ground-state). Eigenstates |0~) constitute a basis where S x operator is diagonal. 
Coming back to the original problem defined by Eqs. (jl2l 1131 ) we note, that in accordance 
with (|T5|) the correct (unique) eigenstate obeying Eq. (fT2j) can be constructed as specific 
linear combination of |0~ ) and |0~ ): 

|G,) = (19) 



which diagonalize operator U\. The eigenstate \Gj) is similar to the eigenstates \G) of 
Ref. Q], cf. Eq.(5) of that paper. 

It follows from the boundary condition (|17|) that shift of the phase 7 by tt does not change 
the boundary problem defined by Eqs. (|16U17|) . Thus the ground-state energy of the system 
and the supercurrent through the circuit are periodic functions of the flux $ c with period 
$o/2. 

At Nw ~ N 2 v / Ej <C 1 fluctuations are weak, the amplitude of potential energy in (|16[) 
is small and its effect is most significant when 4<3> c /$ is integer and the energy levels E m ^ 
are degenerate. In this regime the usual approximation for semiclassical weak link is valid, 
and for the persistent current through the circuit we obtain 

2edE „V2L°7T / 2\~/\\ ( 1 \ . . 

kdl AN V n J \yJ(l-mM* + ( q /2)* J 

Phase- dependent current J (7) is described by equation (|2*0]l for — tt/2 < 7 < n/2 and is a 
periodic function of 7 with period tt. So in the regime of weak fluctuations the dependence 
1(7) demonstrates sawtooth behavior slightly rounded due to fluctuations. 

The opposite limit Nw ^> 1 corresponds to the regime of strong fluctuations. In this case 
the dependence of the eigenvalue b on the phase 7 is exponentially weak : 

b=-2q + 16 W^-g 3 /VV^(l - cos 2j) (21) 
and for the persistent current in the ground state we find 

1(7) = 32 ■ 2 3 / 8 / c ° (v/Ejf 4 VNexp {- 16 - ^ N J^-X sin 27. (22) 



Equation (jzzj), together with Eq.(j9j), presents one of our main results: it gives the amplitude 
of 4e - periodic Josepshson current in the regime of maximal frustration. 

Let us now turn to the investigation of the general situation described by equations (fT3*j) 
and (|12j) . As was mentioned above if the flux per elementary plaquette differs slightly from 
half supercondicting flux quantum the persistent current through the chain has two compo- 
nents l\ e and l2e- In the regime of strong fluctuations both these currents are exponentially 
small. The main exponential factors in the expressions for them can be found on the basis 
of equation (fE^ using the semi-classical approximation. 
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Note that equation ()13j) corresponds to a particle of mass 1 with spin S moving in one- 
dimensional potential 

U(x,S) = wcos2x- S x - hS z , (23) 

the particle energy being E° = b/2. So denoting by 9 and the angels determining the 
spin direction, we can write the imaginary time tunneling amplitude in the form of path 
integral 



(6 l 2 ,02,X 2 |e ™ 1 0i ,02 



d 2 A2,x 2 ( r T/2 f _ ^2 



VilVx exp<- dr[ iS(l - cos 6)<p + — + U(x ,S))} (24) 

{ J-T/2 

For our future purposes it is more convenient to use the above path integral in another form 
also derived in 



5*2 , %2 



e -TH 



Si,x 



C* vSm 5 ( S " - s2 ) ex » { - C dT ( ! ^# + T + ^ ■ *>) } (25) 

We will analyse the expression (f23j) for the limit of relatively large 5 when h ^> w. In 
this regime the field h in (J2*3*j) fixes the direction of S so that S x and S y are always small. 
Therefore we can linearize the action in (}24j) with respect to S x and S y . After linearization 
we can easily exclude the variable S y using the equations of motion. The substitution of 
variables S x — > \fShy, r — > r/h leads to the path integral 

ry2,x2 

(V2 e~ TH \yi ,xi) = / VxVyexp(-S E ) (26) 



where the action 



h r +co 

Se = ~ dr{x 2 + y 2 + U eff (x,y)} , (27) 



U eff (x,y) = (y + dcos2x) 2 + d 2 sm 2 2x , d=^-^- = -^—. (28) 
The appropriate equations of motion are 

x + 2dysin2x = 0, (29) 

y — y — dcos2x = . (30) 
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Using semi-classical approximation we should first determine the classical minima of the 
potential (|2Bj). Within the same limit /i > w we find that U e ff has two groups of minima 
(we call them "even" and "odd" minima) 

x — irn , y = —d , (31) 



7T 

x = - + nn, y = d, (32) 

where n is an arbitrary integer. All these minima correspond to the same value of U e ff = 0. 
So we have to consider two types of tunnelling trajectories. Trajectories of the first type 
connect minima of the same group, i.e. "even-even" and "odd-odd", and corresponding 
variation of the variable x between minima is ±7r, whereas y returns to its original value. 
Trajectories of the second type connect minima of opposite parity (i.e. opposite signs of y), 
and change x variable by ±|. It is not difficult to see from Eqs. 1)111131121 ). that increment 
Ax of the variable x along tunnelling trajectory is in one-to-one correspondence to the 
elementary charge transported along the rhombi chain: go — —Ax. Therefore trajectories of 
the first type lead to 4e - supercurrent, whereas trajectories of the second type produces usual 
2e-quantized supercurrent. The amplitudes of the supercurrent components are determined 
(cf. Appendix for the direct derivation) primarily by the classical actions on corresponding 
trajectories: 

= J 2e sin 7 + I Ae sin(27), (33) 

where 

h e = A 4e exp(-^ e ), J 2e = A 2e exp(-S| e ) , (34) 

and iSjf and Sff are the values of tunnelling actions on trajectories of the first and second 

type respectively. Both S^r and Sjf are large in the region of strong fluctuations Nw 3> 1, 

thus the total supercurrent will in general be dominated by least-action processes. 

To compare actions S E e and we note that the dynamical system (|2l?|) and (|3T)|) has 

two characteristic frequencies. The first one characterises "spin" subsystem with u s — 1, 

whereas the second one is the frequency of the "x" subsystem, uj x ~ d, since typical value 

of y in ()29|) is d. Therefore at d <S 1 - i.e. at sufficiently large flux deflections 5, the "spin 

variable" y is fast and can be integrated out in adiabatic approximation, which leads to 

f ( x 2 d 2 ) 
Se = h / ^ r j"^ J cos 4^ \ , (35) 
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S 4 E e = h I dr{%--d 2 cos2x \ = 4hd (39) 



S% ss 2hd, and S| e « /id, at d < 1 (36) 

The dominant process is thus usual 2e transfer. Comparing the action with the action 
corresponding to the Shrodinger equation (fTB^l 

S° E = J dr^j + q cos2x^j (37) 

and using (J2"T|) we obtain supercurrent amplitude 

/„ « 32 ■ 2-/«tW » ) 3/2 exp (_5^1J^l at J 5 ^«l (38) 

At small flux deflection 5 the parameter d ^> 1 and the spin variable y is relatively slow and 
almost does not change on the type-1 trajectory. The dominant trajectory is then 4e-one. 
Assuming y to be constant, we get 

2 

Taking into account also the first-order term of perturbation theory over 1/d <C 1, we find 
<9|f = /i(4d — 1). Comparison of the action (jHSj) with ()37|1 and ()21|1 allows us to determine 
the pre-exponential factor in the expression for the current 

/^128.2V^(^) 3 ^xp(-^^ + ^} at If »1 (40) 

\Ejy/6j { 7T Ejy/5 7T 2 J 5 3 / 2 £j 

Note that at 5 determined from the equation h = w (where the linear approximation used to 
describe the spin degree of freedom fails), the 4e-current from ()4Uj) matches the exact result 
for 5 = presented in (|2*2*Jl . 

In the intermediate region of d ~ 1 we analyse equations (J27)) . (J29|) and (J3Uj) numerically. 
First we write S E and S"^ e in the form S E = hS E e (d), S E e = hS^i^d). Then the functions 
S E {d) and S E e {d) depending on a single parameter d have been evaluated numerically. 
The result is presented on Fig. El The actions for both types of trajectories are equal 
at d = d f=3 3.2, where we have S'lf(do) = S'ff(do) ~ 11.9. Thus the crossover between 
4e-regime and 2e-regime takes place at 

/ 2 \ 1/3 / \ 2/3 



W = WC =(4^PfJ *° K °- 2 UJ *» (41) 
Varying flux $ r in some vicinity of crossover point (|41|) one can find both 2e and 4e compo- 
nents of supercurrent, but their relative weight is expected to vary strongly with <£> r — 
in some analogy with phase coexistence near first-order phase transition. 
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Figure 3: The results of numerical evaluation of S^r(d) (solid line) and S^{d) (dashed line). 
IV. LOW- VOLTAGE STATES 

In previous Section we obtained estimates (|2(JI22I38I4(J|) for the equilibrium supercurrent 
1(7) around the flux- bias ed rhombi chain with N ^> 1. Note that the maximum value of this 
supercurrent is small, compared to individual critical current of a single junction 7°, even 
in the case of weak quantum fluctuations, cf. Eq.()20|). This is due to the fact that in our 
analysis we have considered perfectly equilibrium Josepshson current, which must be 2n- 
periodic as function of the total phase bias 7. Therefore in the long chain phase differences 
across each rhombi scales as 1/N, leading to I c ~ Ic/N in weak-fluctuation limit Nw <C 1 
(in the opposite limit I c is small exponentially in N). It is clear, however, that under the 
condition of some current bias, with a fixed I -C 1°, the chain will be in some "nearly 
superconducting" state with a very low voltage, due to rare phase slip processes. Below we 
consider regime of relatively large currents (the condition to be specified below), when the 
processes of tunnelling in different rhombi are incoherent. In this case mean voltage V along 
the whole chain can be estimated just as N times the voltage along a single rhombus. Below 
we estimate probability per unit time of an individual QPS in a single rhombus at the fixed 
transport current I <C 1°, and find the V(J) dependence. 

Introducing variables 8 = 8^ + 8^ 2 \ xi = ~ X2 = 8^ 3 ' — 8^ we can rewrite the 
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imaginary-time action for a single rhombus carrying external current / in the form 

1 



St 



dr 



S2E C 



26 2 + xi+X2 2 )+V(8,xi,X2) 



(42) 



V(9, Xi, X2) = -Ej ( 2 cos ^ cos y + 2 sin °- cos y + ^9 



(43) 



We have assumed here that the flux inside the rhombus equals half the sup erco ducting 
flux quantum. In order to find the classical states of the rhombus we eliminate xi an d X2 
from and get 



V{9) 



-E 



.1 



e 



cos ■ 



e 



sin 



(44) 



The potential 
the equation 



has a number of local minima 9 min = 6 + nm where 9 is determined by 

sin y - cos y = — . (45) 

With an appropriate choice of phases xi an d X2 every 9 m i n corresponds to a classical state 
of the rhombus localized near this minimum. Due to quantum tunneling all these states are 
metastable and have finite decay time r. 

Within the semi-classical approximation (valid for Ej ^> Eq) the decay time r is de- 
termined by vicinity of a bounce i.e. a classical trajectory starting at a minimum of the 
potential energy (|43|) coming close to another one and then going back to the first mini- 



mum 

to xi 



111 ]. To be specific we will refer here to the decay rate of a state corresponding 
0, X2 = and 9 = 9q. Decay of this state goes via one of two possible bounce 
trajectories (for I > 0). One of them passes near 9 = 9 , xi — 27r and X2 — while the 
other passes near 9 = 9 , Xi — ~ 2tt and X2 = 0. Both these bounces give equal contribution 
to the decay rate. 

Let us denote by q = (9, xi, X2) T — the three-dimensional column- vector in the coordinate 
space of the rhombus. We also introduce qo(r) = (0q, 0, 0) T as the trajectory corresponding 
to the system being at the minimum of the potential ()43)1 and ^(t) as the bounce trajectory 



which can be determined by solving the classical 
per unit time of the unstable state is given by 



equations of motion. The decay probability 
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Figure 4: Results of numerical evaluation of s(J/iJ). 

where Dei! indicates that the zero eigenvalue is to be omitted when computing the deter- 
minant. 



After changing the time scale according to r — > r/y/EjEc the bounce action can be 



rewritten as S^f^] = 2 \J Ej / Ec s(I/I®) and for the inverse decay time we obtain 



1 n^E c fl^ (TlT0 ,_ .[Ej 



2 V J ^ K(I/I° C ) exp f -2yj^- s(I/I° c ) I . (47) 

Here K(I/I®) is a numerical factor of order one. The function s(I/I®) depending on the only 
parameter I/I® can be evaluated numerically by solving the Lagrangian equations for the 
action (|4~2*|) with the appropriate boundary conditions. The result is presented on Fig. |3J Let 
us assume that the current I is not very small so that the energy difference 5V = ttEjI/I® 
between two nearest minima of the potential (}4~3j) is much larger than the quantum amplitude 
for a phase slip v introduced above, i.e. we assume that I ^> ii = I®v/irEj. In this case 
transitions within each rhombus between the states corresponding to different minima of 
the potential (|43)l are incoherent. Total voltage along the chain can be expressed in terms 
of r as V = Nh6/2e ~ 7cNh/2eT since during each jump of the system from one minimum 
to another the phase 9 changes by n. Thus we obtain for low-current V(I) dependence: 

irNEj 



V(I) 



(^) 1/4 ^(/// c °)ex P (-2y^ S (/// c )) (48) 



15 



0.016 
0.014 
0.012 

O 

o 

c? 0.01 

o 

to 

0.008 

0.006 

0.004 ( 

E / E c 

Figure 5: The critical deviation <5$ c as a function of ratio Ej/Ec 

Equation (jjHj) describes the rhombi chain when the bias current / is large enough: / > I c , 
I ^> I\. Under these conditions coherence in the system is destroyed. This limit is opposite 
to the one we have considered in previous Section, where the value of equilibrium Josephson 
current was determined by coherent quantum fluctuations of all rhombi. 

V. CONCLUSIONS 

In this paper we provide detailed calculations of superconductive current in a long chain 
composed of frustrated rhombi (i.e. loops made of 4 superconductive islands). We show that 
supercurrent carried in 4e quanta dominates over usual 2e supercurrent in the close vicinity 
of the maximally frustrated point $ r = $o/2- According to (}4*T]) the critical deviation <5$ c 
from this point, which brings the system back to usual 2e-supercurrent, depends on the 
only parameter Ej/Ec- This dependence is presented on Fig. We see that 5$ c rapidly 
decreases with the increase of the ratio Ej/Ec- In order to observe experimentally the 4e- 
supprecurent one should control the flux $ r penetrating each rhombus with accuracy better 
than <5$ c . Thus Ej/Ec should not be too large. 

In accordance with ()14j) the parameter q = Nw / 2 governing the strength of fluctuations in 
the maximally frustrated point is proportional to N 2 v/Ej with v defined by equation 0. In 
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Figure 6: The critical current I c at the crossover point as a function of N at differnt Ej/Ec- 

the regime of strong fluctuations (large q) both 4e- and 2e-supercurrents are exponentially 
small, cf. Eq. (J34|) . The actions and .Sjf in ()34|) are proportional to the number of 
rhombi N: Sff ~ N5Sff and Sjf ~ N5S^. At sufficiently large iV small variations of 
and f^S"!; 6 near the crossover point <& c r lead to strong alteration in the relative weight 
of 2e- and 4e-supercurrents. Thus the crossover between 2e- and 4e-regimes is expected to 
be sharp for large N and N 2 v/Ej > 1. On the other hand, the magnitudes of supercurrent 
components J 4e and I 2e , although suppressed by quantum fluctuations, should be not too 
weak to be measured. The semi- qualitative dependence of the critical current I c at the 
crossover point on the number of rhombi at different Ej/Ec is presented on Fig. JBj). While 
calculating the curves depicted on Fig. El the pre-exponential factor in the expression for 
the critical current was evaluated as a geometrical mean of the prefactors in (|38|) and (|4U|) . 
The optimal set of parameters seems to be the following: 5 < Ej/Ec < 7, and iV e (6, 10). 
According to Figs. El and H it would give 5$ c /$ G (0.01 - 0.015) and I c ~ 10" 2 - 10" 3 / c °. 

In this paper we have analysed the path integral (|25j) for the limit of relatively large 8 
when h ^> w. For this condition to hold at the crossover point we need, cf. (|T4| l4T| 



This is always true for large Ej/Eq- However for the proposed set of parameters 
(5 < Ej/E c < 7) the ratio 0.75 < h/w < 0.95 and we are at the edge of the validity 




(49) 
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Figure 7: An analog of dc-SQUID configuration for the measurement of the current-phase rela- 
tion (1331) . 

region for our approximation. Therefore in order to obtain accurate estimates for I 2e and l± e 
at the point of crossover one needs to calculate the classical actions on 2e- and 4e-trajectories 
for the full path integral ([25|) . 

Possible experimental arrangement for testing the current-phase relation (jHSj) is presented 
on Fig. The circuit on Fig. is an analog of a simple dc-SQUID. Let us denote by the 
order parameter phase difference in points A and B. It follows from Eq. (|33j) . that in order 
to evaluate the critical current of the proposed device, much as with a dc-SQUID, one needs 
to maximize over phase the current I given by 

27T$ C \ „ , lT<& r . ( , 7T$ r 



2i 4e cos — ; — sin 



+ 



7T$ ( 

2l2e cos — ; — sm 



(50) 



When the deviation 5$ of the magnetic flux $ r through each rhombus from $o/2 exceeds 
the critical deviation 5$ c , the 4e-supercurrent is negligible and for the critical current of the 
circuit on Fig. (JJJ) we get (in complete analogy with a dc-SQUID) 

7T<£> 
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2e 



COS 



(51) 



So the dependence of 1% on the flux <3> c is $o-P er i°dic for 5$ 3> 5$ c . On the other hand, at 
the maximally frustrated point only 4e-supercurrent survives so that J* is $o/2-per iodic 

2vr$ r 
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•le 



cos • 



$ 



(52) 



The dependence at the crossover point (J 2e = /4 e ) is presented on Fig. [HI 

In our analysis we neglected two intrinsic sources of disorder which are always present 
in the problem considered: a) some weak randomness of fluxes $^ penetrating different 
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Figure 8: The critical current I* at the crossover point. 

rhombi (due to unavoidable differences in their areas), and b) random stray charges q n 
which produce, due to Aharonov-Casher effect, some random phase factors to the phase slip 
tunnelling amplitudes. Whereas the effect of type-a) disorder may be expected to be weak if 
areas of different rhombi coincide with the accuracy better than 5<& c /$o, the b)-type effect 
may occur to be more severe, cf. Ref. [5j, where it was discussed for the simple J J chain. 
We plan to study these effects in further publications. 
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Appendix A: SEMICLASSICAL ANALYSIS FOR I( 7 ) 

In this Appendix we will obtain the dependence of the lowest eigenvalue b of the problem 
defined by Eqs. (ll2ll3l ) on the phase 7 in the regime of strong fluctuations and derive the 
expression (J3*3*j) for the persistent current. 

Let us analyse the transition amplitude (f2T)J) in more details for the case when (xi,yi) is 
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Figure 9: The contour plot of the potential U e ff(x,y). 

an "even" minimum of the potential (}2*Hj) and (22,2/2) is the neares "odd" one, cf. (|31l32j) . 
To be specific we choose (xi,yi) = (0,—d) and (22,2/2) = (7r/2, d). The contour plot of 
the potential U e ff(x,y) is presented on Fig. Possible tunneling trajectories of the system 
are schematically depicted with arrows. It is convenient to divide all trajectories into eight 
groups. Along each trajectory from group 1, variable y is unchanged and equals — d on 
both ends of the trajectory, whereas variable x inreases by it; a trajectory from group 2 is a 
counterpart (going against the arrow on Fig. EJ) to the previous one. The groups 3, 4, . . . , 8 
are defined in the same way according to Fig. |H1 All trajectories from groups 1, 2, 7, 8 
connects minima of the same parity and so are of the first type according to Section 3, 
whereas the trajectories from groups 3, 4, 5, 6 connects minima of opposite parity and are of 
the second type. 

Let us denote by Ta 4e and Ta 2e contributions to the tunneling amplitude from a single 
trajectory of the first and the second type respectively, i.e. 

a 4e = (34 e e~ S ^ and a 2e = P2e e ~ S ^ > (Al) 

where the prefactors /?4 e and /^e can be obtained by integration over the fluctuations near the 
corresponding trajectories. In order to evaluate the transition amplitude (ffijj) in semiclassical 
approximation one should sum up the contributions from all trajectories consisting of n\ 
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subtrajectories from group 1, n 2 subtrajectories from group 2 and so on. Such a trajectory 
including R = Y^k=i Uk subtrajectories gives to the path integral 

rpR 

„ ,n3+n4+n5+n6 _ni+n,2+ri7+n8 / A n\ 

As the trajectory under consideration starts at (0,—d) and ends at (ir/2,d) one should 
impose two additional constraints upon the integers rii, . . . , ng : 

2(n 1 + n 7 ) - 2(n 2 + n 8 ) + n 3 - n 4 + n 5 - n 6 = 1 , (A3) 



n 3 - 77-4 - n 5 + n 6 = 1 . (A4) 

Let us introduce 

K = rii + n-j , L = n 2 + rig and M = n 4 + ri 5 = n 3 + n e — 1 . (A5) 

All trajectories with fixed K, L, M, n 3 and n 4 give the same contribution to the transition 
amplitude. The number of such trajectories is 

^ ! (M + 1)!M! 

(2M + l)\K\L\n 3 W(M - n 4 )!(M - n 3 + 1)! v ; 

Thus, taking into account (jA2IA6IA3|) . for the transition amplitude in semiclassical approx- 
imation we get 

, I -th 1 v \- (Ta ie ) K+L (Ta2e) 2M+1 (M + 1)!M! 

(2/2 ' X2|e l2/1 ' Xl)= 2. K!L!n 3 !n 4 !(2M + l)!(M-n 4 )!(M-n3 + l)! (A?) 

M>max(n4,n3 — 1) 
K — L+n:j—ri4=l 

Instead of calculating the sum in (|A7|) it is convenient to evaluate function Q(a2e,«4e) 
defined by 

, . rp jT , , 1 , . > 

{y2,x 2 \e = ( A8 ) 

For the function Q(a2 e , a4e) we have 

n( s i (r a4e ) K ^(r a2e ) M " 2 i;(Af + 1, m + 1) 

Q( ^ ) = 2 ™ !4! (M-„ 4 )!(M-, !3+ l)! < A9 ' 

M>max (n.4 ,713— 1) 
A'— L+nz— 714=1 

where -B(x, y) is Euler's beta function. Using the integral representation for the beta function 
and the ascending series for the modified Bessel function 7] 

B(x,y)=2j o Msin^-^cos^- 1 , *»(*)= (^*) E fc!( n + A:)! (A1 ° } 
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we can carry out the summation over M, K, L in (|A9|) and obtain 

(Ta 2e )"^+i/ |n4 _ n3+1| (2Ta 4 e 

n3,ri4>0 

Introducing Z = n 3 — 1 — n 4 and accomplishing the summation over n^, n 4 under fixed Z 
we get 



Q= 2^ On 3 +n 4 +l„ |„ j / MB*?) 3+ 4 /|„ 4 -„3+l|(^«2eSm^). (All) 

„ „ 113.114. Jq 



Q = — — V" / rf^/| Z |(2Ta4 e )/|z|(^a2 e sinv9)J| Z+ i|(Ta 2e sinv9) . (A12) 

1 Z=-oc J ° 

Taking into account the integral representation of the modified Bessel function and its gen- 
erating function |j| 

"2n ( _ / 1 \ ^ +°° 



— / <i# exp (m0 + zcostf) , exp || [t + ^ | = ^ t fc 4(z) (A13) 



we can rewrite (jA12|) in the form 

Ta 2e r , f 2n dOxdO: 



Q(a 2e ,a4e) = — r~ / d<P (cos6>i + cos6> 2 )x 

4 Jo Jo W 

exp {Ta 2e sin <^ (cos #1 + cos 6 2 ) + 2Ta 4e cos(6 l i + 9 2 )} . (A14) 

After the substitution u = {61 + 9 2 )/2, v = {61 — 9 2 )/2 we integrate over if and v using the 
relation [3| 

f n / 2 , , , 7r 9 , . sinh 2 z , . * 

/ ^J 1 (2zsin^) = -/ 1 2 /2 (z) = (A15) 

Jo 1 z 

and derive an integral representation for Q(a 2e , a 4e ) 

Q(a 2e ,«4e)= / — {cosh (2Ta 2e cosw) — 1} exp (2Ta 4e cos 2u) . (A16) 
Jo 2tt 

Finally, with the aid of (|A8|) we obtain an explicit expression for the semiclassical transition 
amplitude (J2fi|) 

{y 2 , x 2 \ e~ TH \yi , x\) = / — e m exp(2Ta 4e cos 2u + 2Ta 2e cosu) (A17) 

Jo 2tt 

On the other hand the transition amplitude (J26|) can be written as a sum over the system's 
eigenstates 

(2/2 , x 2 | e~ TH \ Vl ,x 1 ) = J2 yi)^n{x 2 , y 2 )e~ E " T . (A18) 

n 

Note that the symmetry group of the potential U e ft(x,y) consist of transformations 
V n = R n T 7rn / 2} where T a is operator of the translation over distance a along the x-axis 
introduced in Section 3 and R is operator of the reflection in the x-axis. 
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Thus the energy levels E® of a fictitios particle moving in the potential U e ff are classified 
by imposing on their wave functions ip u a twisted boundary condition 

V^ u (x, y) = i) u (x + vr/2, -y) = e iu ip u (x, y) . (A19) 

Comparing (jA18l IA19|) with ()A17|) we conclude that the result (jA17|) has the form of the 
expansion (jA18|) with the multiple e iu under the integral emerging from ip^(xi,yi)ip u (x2,y2) 
and the remaining part of the under-integral expression providing us with the particle energy 

E® = — 2a 4e cos 2u — 2a 2e cos u . (A20) 

Coming back to the original problem defined by Eqs. (fT2l IT3*j) and comparing (fT2*|) 
with (IAT91) we see that we should identify the phase 7 with the " quasimomentum" u. Taking 
into account the relation between b and the energy of the fictitios particle b = 2E° mentioned 
in Section 3 we finally obtain the 6(7) dependence: 

6(7) = — 4ct4 e cos 27 — 4a2e cos 7 . (A21) 

With Eq. (jA21j) and standard relation 1(7) = (2e/h)dE Q /d'y we easily recover the results (JHHj) 
and (H3J). 
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